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Abstract 



We consider a network of cortical spiking neurons with axonal con- 
duction delays and spike timing dependent plasticity (STDP), which is 
ji^ \ representative of a cortical column or hypercolumn. Each neuron fires fol- 

^ . lowing a Hodgkin-Huxley like dynamics and it is interconnected randomly 

to another neurons. In this work we use a permutation entropy approach 
to characterize the dynamic of the spiking neural activity of a given popu- 
lation of neurons. This allow us to investigate the fine temporal structure 
of complex neuronal signals not just by using the probability distributions 
associated with inter-spike intervals, but instead considering much more 
subtle measures accounting for its causal structure. 
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PACS: 02.50.-r , 05.45. Tp , 87.19. La. 
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1 Introduction 

Detection of dynamical changes in complex systems is one of the most relevant 
issues in neuroscience. In the field of experimental neurophysiology, a magni- 
tude used to analyze the network is the normalized variance of the probability 
distribution associated with the inter pikes intervals (ISI). Unfortunately the 
dynamic effects of the network cannot be fully captured by such simple calcula- 
tions, and it does not provide therefore reliable dynamic information that can 
be tested in the phase space. 

One of the simplest computational models of neurons is the Simple Model 
of Spiking Neurons by Izhikevich pLj. This model uses just a two-dimensional 
(2-D) system of ordinary differential equations and four parameters for generate 
twenty different types of neural dynamics. Our research is focused on biophysical 
properties of the Hodgkin- Huxley equations [2 , but using the simple model of 
Izhikevich p] . In order to investigate the dynamical properties of a population 
of neurons, we use a neural network capable of reproducing the main biological 
characteristics of brain oscillation patterns and spontaneous brain activity [3] . 

Brain dynamics variables are noisy, non-stationary, nonlinear and rife with 
temporal discontinuities. Moreover, due to the occurrence of noise and artifacts 
in various forms, it is often not easy to get reliable information from a series of 
measurements. Bandt and Pompe [4] have proposed a robust approach to time 
series analysis on the basis of counting ordinal patterns introducing the concept 
of permutation entropy for quantifying the complexity of a system behind a time 
series. Thus, it is considered the ordinal structure of the time series instead of 
the values themselves [5]. 

In this paper, we use an information-theory-based approach to investigate 
the causality within inter spikes intervals (ISIs). By estimating Fisher infor- 
mation versus MPR complexity [51 [71 |H] of the ISIs signals, we show that it is 
possible to quantify the optimal amount of interconnected neurons that maxi- 
mizes information transmission within a biologically plausible neuronal network 
simulation. This approach provides, thus, dynamic information about the neu- 
ral network that may play a crucial role in the rhythmic activity pattern in the 
brain, and in the processing of information. 

2 Neuronal network simulation model 

Neurons fire spikes when they are near a bifurcation from resting to spiking 
activity, and it is the delicate balance between noise, dynamical currents and 
initial condition what determines the phase diagram of neural activity. While 
there are a huge number of possible ionic mechanisms of excitability and spike 
generation, exist just four bifurcation mechanisms that can result in such a 
transition. That is, there are many ionic mechanisms of spike generation, but 
only four generic bifurcations of equilibria. These bifurcations divide neurons 
into four categories: integrators or resonators, monostable or bistable [9 . 

Bifurcation methodologies [S] allow us to reproduce accurately the biophysi- 



cal properties of Hodgkin- Huxley neuronal models by just taking a two-dimensional 
system of ordinary differential equations and four different parameters. This 
model is named simple model of spiking neurons [I] , and the system of ordinary 
differential equations reads: 

v' = OMv^ + 5v + UO.u + I (1) 

u' ^a{bv~u) (2) 

with the auxiliary after-spike resetting 

if v> +30mV, then J"''- (3) 




Where v is the membrane potential of the neuron, and u is a membrane recovery 
variable, which accounts for the activation of K+ ionic currents and inactivation 
of Na+ ionic currents, and gives negative feedback to v. Thus we are just 
considering a two-dimensional (2-D) system of ordinary differential equations of 
two variables u and v, and all the known types of neurons can be reproduced by 
taking different values of the four parameters a, b, c, d. After the spike reaches 
its apex at -1-30 mV (not to be confused with the firing threshold) the membrane 
voltage and the recovery variable are reset according to Eq. ^. The variable 
/ accounts for the inputs to the neurons [T]. Summarizing, each neuron can 
be described by a simple spiking model that allows us to reproduce the twenty 
most fundamental neuro-computational features of biological neurons [1] . 

Here, we consider a network simulation model in which the number for in- 
terconnected neurons is a parameter under control, and each neuron can be 
interconnected randomly with two neurons or more neurons [3]. The network 
simulation takes into account cortical spiking neurons with axonal conduction 
delays and spike-timing-dependent plasticity (STDP). The magnitude of the 
synaptic weight between pre- and postsynaptic neurons depends on the timing 
of the spikes according to the STDP rules [5^ . Importantly, the interplay be- 
tween conduction delays and STDP helps the spiking neurons to produce stable 
firing patterns that would not be possible without these assumptions. Each 
neuron in the network is described by the simple model of spiking neurons of 
[I] which has been described above in Eq. ^ to Eq. ([3]). 

This model presents biologically plausible phenomena such as patterns of 
spontaneous cortical activity and brain wave like oscillations. Since we can- 
not simulate an infinite-dimensional system on a finite-dimensional lattice, we 
choose a large network with a finite-dimensional approximation taking a time 
resolution of 1 ms. To illustrate the simulation we take first a network consists of 
n — 1000 neurons with n^ = 550 being of excitatory regular (RS) spiking type, 
and the remaining n^ — 450 of inhibitory fast spiking (FS) type {n = n^ + Ui). 
Each excitatory regular spiking type neuron is connected to m — 2,6,8, 10, 20, 
30, 40, 50, 60, 70, 80, 90, 100, and 120 random neurons, so that the probability 
of connection is m/n = 0.002, 0.006, 0.008,0.01,0.02,0.04,0.06,0.08,0.1 and 0.12. 
Each inhibitory fast (FS) spiking neuron is connected to to = 2,6,8, 10, 20, 



30, 40, 50, 60, 70, 80, 90, 100, and 120 excitatory neurons only. The second a 
network consists of n = 900 neurons with A^e = 500 being of excitatory (RS) 
type, and the remaining ni = 400 of inhibitory (FS) type (n = rig + n^). In 
this case, each excitatory regular spiking type neuron is connected to to = 2, 
4, 6, 8, 10, 20,30,40,60,80,100 and 120 random neurons, so that the probabiUty 
of connection is m/n = 0.002, 0.004,0.006,0.008,0.01,0.02,0.03,0.04,0.06,0.08,0.1 
and 0.12. And each inhibitory (FS) neuron is connected to r?T- = 2, 4, 6, 8, 10, 
20,30,40,60,80,100 and 120 excitatory neurons only. In both cases synaptic con- 
nections among neurons have fixed conduction delays, which are random integers 
between 1 ms and 10 ms. The two network arrangement described above are 
representative of the general behavior shown by large network simulations tak- 
ing the axonal conduction delays and spike-timing-dependent plasticity (STDP) 
approach. Although the pattern are random, reflecting connectivity within the 
cortex, one can implement sophisticated anatomy and the two sparse network 
being shown in this section are representative of a cortical column or hypercol- 
umn [5]. 

3 Information Theory quantifiers 

3.1 Shannon entropy, Fisher Information Measure and MPR 
- Statistical Complexity 

Sequences of measurements (or observations) constitute the basic elements for 
the study of natural phenomena. In particular, from these sequences, com- 
monly called time series, one should judiciously extract information on dynami- 
cal systems under study. We can defined an Information Theory quantifier as a 
measure that is able to characterize some property of the probability distribu- 
tion function (PDF) associated with these time series. Entropy, regarded as a 
measure of uncertainty, is the most paradigmatic example of these quantifiers. 
Given a time series X{t) = {xt; t = I,--- ,M}, a set of AI measures of the 
observable X and the associated PDF, given by P = {pj;j — 1, • • • ,N} with 
^7=1 Pj = 1 9,nd N the number of possible states of the system under study, 
the Shannon's logarithmic information measure jlO| is defined by 

N 

S[P] = -Y,P,Hp,). (4) 

This functional is equal to zero when we are able to predict with full certainty 
which of the possible outcomes j whose probabilities are given by pj will actually 
take place. Our knowledge of the underlying process, described by the proba- 
bility distribution, is maximal in this instance. In contrast, this knowledge is 
minimal, commonly, for a uniform distribution P^ — {pj — 1/N, Vj = 1, • • • , N}. 
The Shannon entropy S* is a measure of "global character" that it is not too 
sensitive to strong changes in the PDF taking place on small-size region. Such 
is not the case with the Fisher's Information Measure F[f] — J {|V/p//} dx 



[TTl [T2] , which constitutes a measure of the gradient content of the distribution 
/ (continuous PDF), thus being quite sensitive even to tiny locahzed pertur- 
bations. The Fisher's Information Measure can be variously interpreted as a 
measure of the abihty to estimate a parameter, as the amount of information 
that can be extracted from a set of measurements, and also as a measure of 
the state of disorder of a system or phenomenon [T^l [13] , being its most im- 
portant property the so-called Cramer-Rao bound. It is important to remark 
that the gradient operator significantly influences the contribution of minute 
local /-variations to Fisher information value, so that the quantifier is called a 
"local" one. Note that Shannon's entropy decreases with skewed distribution, 
while Fisher's information increases in such a case. Local sensitivity is useful in 
scenarios whose description necessitates appeal to a notion of "order" O [71 [8] . 
The concomitant problem of loss of information due to the discretization has 
been thoroughly studied (see, for instance, [TH [THl [H] and references therein) 
and, in particular, it entails the loss of Fisher's shift-invariance, which is of no 
importance for our present purposes. 

For the Fisher Information Measure computation (discrete PDF) we follow 
the proposal of Ferri and coworkers [17] (among others) 

F[P] = ] V 2^^^±1^ . (5) 

If our system is in a very ordered state and thus is represented by a very narrow 
PDF, we have a Shannon Entropy S" ~ and a Fisher's Information Measure 
F ^ Ffnax- On the other hand, when the system under study lies in a very 
disordered state one gets an almost flat PDF and S ^ Smax while F ~ 0. Of 
course, Smax and Fmax are, respectively, the maximum values for the Shannon 
Entropy and Fisher Information Measure. One can state that the general be- 
havior of the Fisher Information Measure is opposite to that of the Shannon 
Entropy [18]. 

It is well known, however, that the degree of structure present in a pro- 
cess is not quantifled by randomness measures and, consequently, measures of 
statistical or structural complexity are necessary for a better understanding 
(characterization) the system dynamics represented by their time series [19] . 
The opposite extremes of perfect order (i.e.: a periodic sequence) and maximal 
randomness (i.e.: a fair coin toss) are very simple to describe because they do 
not have any structure. The complexity should be zero in these cases. At a 
given distance from these extremes, a wide range of possible degrees of physi- 
cal structure exists. The complexity measure allows one to quantify this array 
of behavior [20] . In the present work we have considered the effective statisti- 
cal complexity measure (SCM) introduced by Lamberti et al. [21], (also called 
MPR-complexity) since it is able to detect essential details of the dynamics. 

This statistical complexity measure is deflned, following the seminal notion 
advanced by Lopez-Ruiz et al. [22| . through the product 

c.,s[P] = Q.j[P,Pe]-ns[P] (6) 



of the normalized Shannon entropy 

Hs[P] = S[P]/S„,a. (7) 

with Smax — S[Pe] — IniV, (0 < Hs < 1) and the disequihbrium Qj defined in 
terms of the Jensen-Shannon divergence. That is, 

Qj[P,Pe] =QoJ[P,Pe] (8) 

with 

J[PPe] = S[{P + Pe)/2]-S[P]/2-S[Pe]/2 (9) 

the above-mentioned Jensen-Shannon divergence and Qq a normalization con- 
stant (0 < Qj < 1), equal to the inverse of the maximum possible value of 
J^[P,Pe]. This value is obtained when one of the components of P, say Pm, 
is equal to one and the remaining pj are equal to zero. The Jensen-Shannon 
divergence, that quantifies the difference between two (or more) probability dis- 
tributions, is especially useful to compare the symbolic composition between 
different sequences [53]. Note that the above introduced SCM depends on two 
different probability distributions, the one associated to the system under anal- 
ysis, P, and the uniform distribution, P^. Furthermore, it was shown that for a 
given value of "Hg, the range of possible Cjs values varies between a minimum 
Cmin and a maximum Cmax, restricting the possible values of the SCM in a given 
complexity-entropy plane |24j . Thus, it is clear that important additional in- 
formation related to the correlational structure between the components of the 
physical system is provided by evaluating the statistical complexity measure. 

3.2 The Bandt-Pompe approach to the PDF determina- 
tion 

The study and characterization of time series X{t) by recourse to Information 
Theory tools assume that the underlying PDF is a priori given. Per contra, part 
of the concomitant analysis involves extracting the PDF from the data and there 
is no univocal procedure with which everyone agrees. Bandt and Pompe (BP), 
almost ten years ago, introduced a successful methodology for the evaluation of 
the PDF associated to scalar time series data using a symbolization technique 
[3]. For a didactic description of the approach, as well as, its main biomedical 
and econophysics applications, see [51. 

The pertinent symbolic data are (i) created by ranking the values of the 
series and (ii) defined by reordering the embedded data in ascending order, 
which is tantamount to a phase space reconstruction with embedding dimen- 
sion (pattern-length) D and time lag r. In this way it is possible to quantify the 
diversity of the ordering symbols (patterns) derived from a scalar time series. 
Note that the appropriate symbol-sequence arises naturally from the time series 
and no model-based assumptions are needed. In fact, the necessary "partitions" 
are devised by comparing the order of neighboring relative values rather than 
by apportioning amplitudes according to different levels. This technique, as 



opposed to most of those in current practice, takes into account the temporal 
structure of the time series generated by the physical process under study. This 
feature allows us to uncover important details concerning the ordinal structure 
of the time series [351 IHl [2S] and can also yield information about temporal cor- 
relation [571 HH] • It is clear that this type of analysis of time series entails losing 
some details of the original series' amplitude-information. Nevertheless, by just 
referring to the series' intrinsic structure, a meaningful difhculty reduction is in- 
deed achieved by Bandt and Pompe with regards to the description of complex 
systems. The symbolic representation of time series by recourse to a comparison 
of consecutive (r — 1) or non-consecutive (r > 1) values allows for an accurate 
empirical reconstruction of the underlying phase-space, even in the presence 
of weak (observational and dynamical) noise [4]. Furthermore, the ordinal- 
pattern's associated PDF is invariant with respect to nonlinear monotonous 
transformations. Accordingly, nonlinear drifts or scalings artificially introduced 
by a measurement device will not modify the quantifiers' estimation, a nice 
property if one deals with experimental data (see, i.e., [IH])- These advantages 
make the Bandt and Pompe methodology more convenient than conventional 
methods based on range partitioning (i.e.: PDF based on histogram). 

Additional advantages of the method reside in (i) its simplicity (we need few 
parameters: the pattern length/embedding dimension D and the embedding de- 
lay r and (ii) the extremely fast nature of the pertinent calculation- process [3Dj . 
The BP methodology can be applied not only to time series representative of 
low dimensional dynamical systems but also to any type of time series (regular, 
chaotic, noisy, or reality based). In fact, the existence of an attractor in the 
D-dimensional phase space is not assumed. The only condition for the applica- 
bility of the Bandt-Pompe methodology is a very weak stationary assumption 
(that is, for k < D, the probability for xt < Xt+k should not depend on t [1]). 

To use the Bandt and Pompe [4] methodology for evaluating the PDF P 
associated to the time series (dynamical system) under study, one starts by 
considering partitions of the pertinent D-dimensional space that will hopefully 
"reveal" relevant details of the ordinal-structure of a given one-dimensional time 
series X{t) = {xt] i = 1, • • • , M} with embedding dimension I? > 1 (Z? G N) and 
embedding time delay r (r S N). We are interested in "ordinal patterns" of order 
(length) D generated by (s) ^ {xs-(D-i)r, Xs-(D-2)r, ■■■ , Xs-r, Xs ), 
which assigns to each time s the ZJ-dimensional vector of values at times s, s — 
r, • • • ,s — {D — 1)t. Clearly, the greater the _D— value, the more information 
on the past is incorporated into our vectors. By "ordinal pattern" related to 
the time (s) we mean the permutation tt = (rp, ri, ■ ■ ■ , rjj^i) of [0, 1, ■ • ■ , D — 
1] de&ned by Xs-ro-iT < Xs^ru-2T < ■■■ < Xs-nr < Xs^ror- In 
order to get a unique result we set r^ < ri_i if Xg-n = Xg-ri-i- This is 
justified if the values of Xt have a continuous distribution so that equal values 
are very unusual. Thus, for all the D\ possible permutations tt of order D, 
their associated relative frequencies can be naturally computed by the number 
of times this particular order sequence is found in the time series divided by the 
total number of sequences. 

Consequently, it is possible to quantify the diversity of the ordering sym- 



bols (patterns of length D) derived from a scalar time-series, by evaluating 
the so-called permutation entropy, the permutation statistical complexity and 
permutation Fisher information measure (that Shannon's entropy, MPR - Sta- 
tistical entropy and Fisher information measure special versions corresponding 
to the Bandt and Pompe-PDF). Of course, the embedding dimension D plays 
an important role in the evaluation of the appropriate probability distribution 
because D determines the number of accessible states D\ and also conditions 
the minimum acceptable length M ^ D! of the time series that one needs in 
order to work with a reliable statistics [5T| . 

With respect to the selection of the other parameter, Bandt and Pompe 
suggest to work with 4 < D < 6 and specifically considered an embedding delay 
r = 1 in their cornerstone paper [4]. Nevertheless, it is clear that other values 
of r could provide additional information. It has been recently shown that this 
parameter is strongly related, if it is relevant, with the intrinsic time scales of 
the system under analysis [311 1331 131] ■ 

3.3 Causal information planes 

The above introduced based Information Theory quantifiers Us, J^ and Cjs 
evaluated using Bandt and Pompe's PDF allow us to define three causality in- 
formation planes: HxC, Hx J^ and C x J^. The causality Shannon-complexity 
plane, H x C, is based only on global characteristics of the associated time 
series Bandt and Pompe PDF (both quantities are defined in terms of Shan- 
non entropies), while the causality Shannon-Fisher plane, T-L x T and causality 
complexity-Fisher, C x F , are based on global and local characteristics of the 
PDF. In the case oiT-LxC the variation range is [0, 1] x [CmimCmax] (with Cmm 
and Cmax the minimum and maximum statistical complexity values, respectively, 
for a given T-Ls value [M]), while in the causality H x T and C x T the range is 
[0, 1] X [0, 1] in both cases. These causal information planes have been profitable 
used to separate and differentiate amongst chaotic and deterministic systems 
|25l [8] ; visualization and characterization of different dynamical regimes when 
the systems parameters vary [6l[7ll3; time dynamical evolution 31,; identifying 
periodicities in natural time series [35] ; identification of deterministic dynamics 
contaminated with noise [361137) and; estimating intrinsic time scales of delayed 
systems |32j |331 [31] ; among other applications (see [5i and references therein) . 

4 Results and Discussion 

In order to illustrate our method we run the simulation for 20000 ms taking a 
time resolution of 1 ms, which is enough to guarantee the condition M >> Dl to 
be satisfied. Figures [T] A and B show the spike rasters considering N — 1000 and 
N = 900, respectively, for a time window of 1000ms. Figures [2] A and B show 
the interspikes intervals (ISIs) considering N — 1000 and N = 900, respectively, 
for the entire time we run the simulation. Note that the previous plots only 
confirm that the neurons within the model fires with millisecond precision and 



we cannot gain any significative information about the dynamic. Figures [3] A 
and B show the normahzed variance Rp = ^JVar{tp)/ < ip > of the ISIs for two 
different simulated systems. Note that increasing degree of inter-connectivity m 
imphes greater data dispersion. Thus, this kind of analysis does not provide us 
any insight of how information might saturates as the degree of interconnectivy 
becomes bigger. 

In the following we use the Bandt and Pompe [4] methodology for evaluating 
the PDF P associated to the time series, considering an embedded dimension- 
ality £) = 6 (with T = 1). This embedded dimensionality is enough to capture 
efficiently causality information of the ordinal structure of the time series [3]. 
Figures [3] A and B show Fisher information Eq.® versus the MPR complexity 
Eq.® for the same biophysical parameters used above. Notice that Figures |4] A 
and B present a maximum at ?7i = 50 and m = 30, respectively, for n — 1000 and 
n = 900. A similar behavior can be observed in Figures [S] A and B, which show 
Fisher information Eq.([S|) versus permutation entropy Eq.Q. Consequently, 
building a phase space Fisher versus MPR complexity (or Fisher versus permu- 
tation entropy) provides us an useful tool to detect and quantify the optimal 
biophysical values which allows the neural network to transmit information more 
efficiently. 

A generic neuron receives inputs from over ten thousand other neurons 
through dendritic connections called synapses. These inputs produce trans- 
membrane currents changing the membrane potential of the neuron. Synaptic 
currents are called postsynaptic potentials (PSP). Small currents produce small 
PSP; large currents cause significant PSP to be amplified by sensitive channels, 
and may change the voltage in the cell membrane. This can leads therefore to 
the generation of a action potential. While there is a great amount of ionic 
mechanisms that may explain the excitability and generation of firing an spike, 
there are only four bifurcation mechanisms that may explain the transition. If 
we consider the geometry of the phase space in these branches, we can under- 
stand many properties computational neurons, and the nature of the threshold. 
We can understand also how these properties relate to each other, why some 
are equivalent and mutually exclusive (see [5] for an exhaustive review). 

Here, we use an alternative methodology when considering a large popula- 
tion of neurons, we shown that by building the phase space of Fisher information 
versus MPR complexity /permutation entropy it is possible to characterize the 
dynamic of the system, and to find the optimal values which maximize infor- 
mation transmission within a neuronal network. 

It is expected that information should saturate when the number of neurons 
increases [38l|39]. However, very little is known about how it could saturate. Our 
current approach provides a novel methodology to respond how deal saturation 
with information taking a realistic neuronal network. The approach is not only 
important from a theoretical point of view; it might help determine which areas 
of the cortex could have a higher level of information, and to evaluate how causal 
interactions in neural dynamics would be modulated by behavior. We believe 
that this will become a very important tool for future research on the encoding 
capacity of biologically realistic neural networks. 
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Figure 1: Spike rasters. A) n — 
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Figure 2: Histograms of the interspike intervals (ISIs). A) ISI for n = 1000, 
He = 550, n, = 450 and m = 50. B) ISI for n = 900, Ue = 500, n, = 400 and 
m = 30. 
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Figure 3: Normalized variance of the ISIs. A) n — 1000, iie ~ 550, Ui — 450, 
and the points correspond to the different interconnectivities m = 2, 6, 8, 10, 
20, 30, 40, 50, 60, 70, 80, 90, 100, and 120. B) n = 900, n^ = 500, m = 400, 
and the points correspond to the different interconnectivities ?ti = 2, 4, 6, 8, 10, 
20,30,40,60,80,100 and 120. 
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Figure 4: Fisher information versus MPR complexity, embedded dimensionality 
D = 6. A) n — 1000, Ue = 550, Ui — 450 and the points correspond to the 
different interconnectivities m = 2,6,8, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 
and 120 (m gets bigger as the variable on the X-axis increases). The maximum 
value of Fisher information corresponds to to = 50. B) n = 900, Ue = 500, 
rii = 400 and the points correspond to the different interconnectivities m = 2, 
4, 6, 8, 10, 20,30,40,60,80,100 and 120 (to gets bigger as the variable on the 
X-axis increases). The maximum value of Fisher information corresponds to 
m = 30. 
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Figure 5: Fisher information versus normalized Shannon permutation entropy, 
embedded dimensionahty D = 6. A) n = 1000, n^ ~ 550, n^ — 450 and the 
maximum value of Fisher information corresponds to n = 50. Same curve labels 
as in Fig. |4]A. B) n — 900, rig = 500, Ui — 400 and the maximum value of 
Fisher information corresponds to to = 30. Same curve labels as in Figures |4]B. 
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